options(scipen = 9999)

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
set.seed(415)

library(rstan)
library(shinystan)
library(stringr)
library(dplyr)
library(readxl)
library(gtools)
library(rjags)
library(runjags)
library(brms)
library(knitr)
library(DT)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

R Markdown

candidates <- c("PS", "CHEGA", "PCP", "PSD", "BE", "IL", "RIR")


pt_pitagorica <- read_excel("data/pollsportugal.xltx", sheet = "Pitagorica")%>% 
  filter(Party %in% candidates)%>% 
  mutate(Bias = Result  - Poll ) %>% 
  mutate(Company = "Pitagorica")


pt_UCP <- read_excel("data/pollsportugal.xltx", sheet = "UCP")%>% 
  filter(Party %in% candidates)%>% 
  mutate(Bias = Result  - Poll ) %>% 
  mutate(Company = "UCP")


pt_eurosondagem <- read_excel("data/pollsportugal.xltx", sheet = "Eurosondagem")%>% 
  filter(Party %in% candidates)%>% 
  mutate(Bias = Result  - Poll ) %>% 
  mutate(Company = "Eurosondagem")


pt_intercampus <- read_excel("data/pollsportugal.xltx", sheet = "Intercampus")%>% 
  filter(Party %in% candidates)%>% 
  mutate(Bias = Result  - Poll ) %>% 
  mutate(Company = "Intercampus")


pt_aximage <- read_excel("data/pollsportugal.xltx", sheet = "Aximage")%>% 
  filter(Party %in% candidates)%>% 
  mutate(Bias = Result  - Poll ) %>% 
  mutate(Company = "Aximage")



pt_bias_final <- rbind(pt_pitagorica, pt_UCP, pt_eurosondagem, pt_intercampus, pt_aximage) %>% 
  mutate(Companyiid = match(Company , sort(unique(Company))))  %>% 
  mutate(Eleciid = match(Election , sort(unique(Election)))) %>% 
  mutate(Partyiid = match(Party , sort(unique(Party))))

DT::datatable(pt_bias_final)
boxplot(Bias ~ Company, data=pt_bias_final, main="Boxplot Bias all elections portuguese companies")

boxplot(Bias ~ Company, data=pt_bias_final %>% filter(Election == "Pres"), main="Boxplot Bias President portuguese companies")

Model with just house and party

mod_string = " 
model {
## sampling
for (i in 1:N){
   y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] , invsigma2)
}

## priors
for (j in 1:J){
   mu_party[j] ~ dnorm(mu, invtau2)
}

for (h in 1:H){
   mu_house[h] ~ dnorm(mu, invtau2)
}


invsigma2 ~ dgamma(1.5, 6)
sigma <- sqrt(pow(invsigma2, -1))

## hyperpriors

mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))

}"





params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")



the_data <- list("y" = pt_bias_final$Bias,
                 "PartyIndex" = pt_bias_final$Partyiid,
                 "HouseIndex" = pt_bias_final$Companyiid,
                 "N" = length(pt_bias_final$Bias),
                 "J" = length(unique(pt_bias_final$Partyiid)),
                 "H" = length(unique(pt_bias_final$Companyiid))
                 )
params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")


mod = jags.model(textConnection(mod_string), data=the_data, n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 69
##    Unobserved stochastic nodes: 14
##    Total graph size: 260
## 
## Initializing model
update(mod, 1e3)

mod_sim = coda.samples(model=mod,
                       variable.names=params,
                       n.iter=5e3)

mod_csim = as.mcmc(do.call(rbind, mod_sim))

plot(mod_sim)

traceplot(mod_sim, col = 1:3)

gelman.diag(mod_sim)
## Potential scale reduction factors:
## 
##             Point est. Upper C.I.
## invsigma2            1          1
## invtau2              1          1
## mu                   1          1
## mu_house[1]          1          1
## mu_house[2]          1          1
## mu_house[3]          1          1
## mu_house[4]          1          1
## mu_house[5]          1          1
## mu_party[1]          1          1
## mu_party[2]          1          1
## mu_party[3]          1          1
## mu_party[4]          1          1
## mu_party[5]          1          1
## mu_party[6]          1          1
## 
## Multivariate psrf
## 
## 1
dic = dic.samples(mod, n.iter=1e3)

(pm_params = colMeans(mod_csim))
##   invsigma2     invtau2          mu mu_house[1] mu_house[2] mu_house[3] 
##  0.18529063  2.83215636 -0.28701191 -0.22371688 -0.34520823 -0.29310483 
## mu_house[4] mu_house[5] mu_party[1] mu_party[2] mu_party[3] mu_party[4] 
##  0.09527954 -0.91251443  0.16557152 -0.18953359 -0.12248021  0.05642923 
## mu_party[5] mu_party[6] 
## -1.07897684 -0.79943599
summary(mod_sim)
## 
## Iterations = 1001:6000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                 Mean     SD  Naive SE Time-series SE
## invsigma2    0.18529 0.0325 0.0002654       0.000294
## invtau2      2.83216 1.4721 0.0120197       0.021036
## mu          -0.28701 0.2411 0.0019683       0.003045
## mu_house[1] -0.22372 0.5045 0.0041192       0.004898
## mu_house[2] -0.34521 0.4972 0.0040597       0.004686
## mu_house[3] -0.29310 0.5229 0.0042698       0.004971
## mu_house[4]  0.09528 0.5499 0.0044898       0.005681
## mu_house[5] -0.91251 0.4113 0.0033586       0.004418
## mu_party[1]  0.16557 0.4729 0.0038609       0.005071
## mu_party[2] -0.18953 0.6314 0.0051553       0.005968
## mu_party[3] -0.12248 0.6261 0.0051121       0.005772
## mu_party[4]  0.05643 0.6260 0.0051116       0.006400
## mu_party[5] -1.07898 0.4620 0.0037724       0.004992
## mu_party[6] -0.79944 0.4364 0.0035635       0.004503
## 
## 2. Quantiles for each variable:
## 
##                2.5%     25%      50%      75%    97.5%
## invsigma2    0.1281  0.1623  0.18338  0.20583  0.25516
## invtau2      0.8791  1.7814  2.53183  3.55290  6.48869
## mu          -0.7461 -0.4490 -0.29317 -0.12875  0.20688
## mu_house[1] -1.1972 -0.5587 -0.23010  0.10560  0.79556
## mu_house[2] -1.3202 -0.6709 -0.35182 -0.02563  0.66141
## mu_house[3] -1.3200 -0.6375 -0.29779  0.04362  0.76527
## mu_house[4] -0.9287 -0.2806  0.07539  0.44576  1.22297
## mu_house[5] -1.7501 -1.1803 -0.89741 -0.63392 -0.12991
## mu_party[1] -0.7158 -0.1585  0.15043  0.47293  1.12995
## mu_party[2] -1.4283 -0.6053 -0.20442  0.20857  1.09561
## mu_party[3] -1.3096 -0.5367 -0.14784  0.27320  1.18168
## mu_party[4] -1.0920 -0.3611  0.02871  0.43778  1.37623
## mu_party[5] -2.0235 -1.3776 -1.06773 -0.76618 -0.20058
## mu_party[6] -1.6898 -1.0789 -0.78776 -0.50435  0.02468
aaa <- summary(mod_sim)
# xtable(as.data.frame(aaa$statistics))

Model with house, party and election type

mod_string2 = " 
model {
## sampling
for (i in 1:N){
   y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] + mu_elec[ElecIndex[i]] , invsigma2)
}

## priors
for (j in 1:J){
   mu_party[j] ~ dnorm(mu, invtau2)
}

for (h in 1:H){
   mu_house[h] ~ dnorm(mu, invtau2)
}

for (r in 1:R){
   mu_elec[r] ~ dnorm(mu, invtau2)
}


invsigma2 ~ dgamma(1.5, 6)
sigma <- sqrt(pow(invsigma2, -1))

## hyperpriors

mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))

}"



set.seed(115)


params = c("mu_party", "mu_house","mu_elec", "mu", "invtau2", "invsigma2")

the_data2 <- list("y" = pt_bias_final$Bias,
                 "PartyIndex" = pt_bias_final$Partyiid,
                 "HouseIndex" = pt_bias_final$Companyiid,
                 "ElecIndex" = pt_bias_final$Eleciid,
                 "N" = length(pt_bias_final$Bias),
                 "J" = length(unique(pt_bias_final$Partyiid)),
                 "H" = length(unique(pt_bias_final$Companyiid)),
                 "R" = length(unique(pt_bias_final$Eleciid))
)
mod2 = jags.model(textConnection(mod_string2), data=the_data2, n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 69
##    Unobserved stochastic nodes: 18
##    Total graph size: 354
## 
## Initializing model
update(mod2, 1e3)

mod_sim2 = coda.samples(model=mod2,
                       variable.names=params,
                       n.iter=5e3)

mod_csim2 = as.mcmc(do.call(rbind, mod_sim2))

traceplot(mod_sim2, col = 1:3)

gelman.diag(mod_sim2)
## Potential scale reduction factors:
## 
##             Point est. Upper C.I.
## invsigma2            1          1
## invtau2              1          1
## mu                   1          1
## mu_elec[1]           1          1
## mu_elec[2]           1          1
## mu_elec[3]           1          1
## mu_elec[4]           1          1
## mu_house[1]          1          1
## mu_house[2]          1          1
## mu_house[3]          1          1
## mu_house[4]          1          1
## mu_house[5]          1          1
## mu_party[1]          1          1
## mu_party[2]          1          1
## mu_party[3]          1          1
## mu_party[4]          1          1
## mu_party[5]          1          1
## mu_party[6]          1          1
## 
## Multivariate psrf
## 
## 1
dic = dic.samples(mod2, n.iter=1e3)


(pm_params = colMeans(mod_csim2))
##   invsigma2     invtau2          mu  mu_elec[1]  mu_elec[2]  mu_elec[3] 
##  0.22210356  2.06050330 -0.30899280 -1.45596782 -0.25251530  0.67944498 
##  mu_elec[4] mu_house[1] mu_house[2] mu_house[3] mu_house[4] mu_house[5] 
## -0.46195493 -0.22388919 -0.42778846 -0.25185377  0.06816059 -0.99148698 
## mu_party[1] mu_party[2] mu_party[3] mu_party[4] mu_party[5] mu_party[6] 
##  0.20019166 -0.35588102 -0.26834499  0.35541802 -1.25080284 -0.80848883
summary(mod_sim2)
## 
## Iterations = 1001:6000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                 Mean      SD  Naive SE Time-series SE
## invsigma2    0.22210 0.04162 0.0003399      0.0004504
## invtau2      2.06050 1.13126 0.0092367      0.0203188
## mu          -0.30899 0.21647 0.0017675      0.0023309
## mu_elec[1]  -1.45597 0.79879 0.0065221      0.0104836
## mu_elec[2]  -0.25252 0.58265 0.0047573      0.0066304
## mu_elec[3]   0.67944 0.49682 0.0040565      0.0083672
## mu_elec[4]  -0.46195 0.48444 0.0039554      0.0066320
## mu_house[1] -0.22389 0.54322 0.0044354      0.0060911
## mu_house[2] -0.42779 0.54539 0.0044531      0.0065991
## mu_house[3] -0.25185 0.56918 0.0046473      0.0061467
## mu_house[4]  0.06816 0.58492 0.0047759      0.0060925
## mu_house[5] -0.99149 0.46061 0.0037609      0.0065718
## mu_party[1]  0.20019 0.49779 0.0040644      0.0061182
## mu_party[2] -0.35588 0.68940 0.0056289      0.0064711
## mu_party[3] -0.26834 0.67890 0.0055432      0.0062635
## mu_party[4]  0.35542 0.70350 0.0057440      0.0081299
## mu_party[5] -1.25080 0.52259 0.0042669      0.0077507
## mu_party[6] -0.80849 0.47878 0.0039092      0.0063874
## 
## 2. Quantiles for each variable:
## 
##                2.5%     25%     50%      75%   97.5%
## invsigma2    0.1486  0.1928  0.2191  0.24887  0.3109
## invtau2      0.6264  1.2565  1.8109  2.59124  4.9546
## mu          -0.7307 -0.4505 -0.3108 -0.16725  0.1237
## mu_elec[1]  -3.1928 -1.9394 -1.3959 -0.91002 -0.0566
## mu_elec[2]  -1.3786 -0.6383 -0.2645  0.12766  0.9244
## mu_elec[3]  -0.2374  0.3409  0.6546  0.99734  1.7289
## mu_elec[4]  -1.4319 -0.7807 -0.4636 -0.15055  0.5043
## mu_house[1] -1.3102 -0.5769 -0.2185  0.12886  0.8571
## mu_house[2] -1.5252 -0.7821 -0.4180 -0.06726  0.6295
## mu_house[3] -1.3764 -0.6260 -0.2546  0.12432  0.8514
## mu_house[4] -1.0787 -0.3149  0.0578  0.45013  1.2319
## mu_house[5] -1.9251 -1.2882 -0.9765 -0.68052 -0.1219
## mu_party[1] -0.7609 -0.1285  0.1923  0.52150  1.2002
## mu_party[2] -1.7513 -0.7901 -0.3474  0.08937  0.9751
## mu_party[3] -1.6207 -0.7015 -0.2706  0.16700  1.0873
## mu_party[4] -0.9384 -0.1253  0.3136  0.80061  1.8541
## mu_party[5] -2.3336 -1.5840 -1.2307 -0.89727 -0.2795
## mu_party[6] -1.7774 -1.1211 -0.7978 -0.48500  0.1082
bbb <-summary(mod_sim2)

# xtable(as.data.frame(bbb$statistics))

Sensitivity analysis

For the sensitivity analysis we will assume very high values for our uncertain variance and nu_o

Model with just house and party

mod_string = " 
model {
## sampling
for (i in 1:N){
   y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] , invsigma2)
}

## priors
for (j in 1:J){
   mu_party[j] ~ dnorm(mu, invtau2)
}

for (h in 1:H){
   mu_house[h] ~ dnorm(mu, invtau2)
}


invsigma2 ~ dgamma(1.5, 15)
sigma <- sqrt(pow(invsigma2, -1))

## hyperpriors

mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))

}"





params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")



the_data <- list("y" = pt_bias_final$Bias,
                 "PartyIndex" = pt_bias_final$Partyiid,
                 "HouseIndex" = pt_bias_final$Companyiid,
                 "N" = length(pt_bias_final$Bias),
                 "J" = length(unique(pt_bias_final$Partyiid)),
                 "H" = length(unique(pt_bias_final$Companyiid))
                 )


params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")


mod = jags.model(textConnection(mod_string), data=the_data, n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 69
##    Unobserved stochastic nodes: 14
##    Total graph size: 260
## 
## Initializing model
update(mod, 1e3)

mod_sim = coda.samples(model=mod,
                       variable.names=params,
                       n.iter=5e3)

mod_csim = as.mcmc(do.call(rbind, mod_sim))


(pm_params = colMeans(mod_csim))
##   invsigma2     invtau2          mu mu_house[1] mu_house[2] mu_house[3] 
##  0.17612451  2.86514613 -0.29513064 -0.23349843 -0.36775627 -0.29418506 
## mu_house[4] mu_house[5] mu_party[1] mu_party[2] mu_party[3] mu_party[4] 
##  0.06359648 -0.90546652  0.15046624 -0.19755661 -0.13760190  0.03354543 
## mu_party[5] mu_party[6] 
## -1.05207361 -0.78294849
mod_string = " 
model {
## sampling
for (i in 1:N){
   y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] , invsigma2)
}

## priors
for (j in 1:J){
   mu_party[j] ~ dnorm(mu, invtau2)
}

for (h in 1:H){
   mu_house[h] ~ dnorm(mu, invtau2)
}


invsigma2 ~ dgamma(1.5, 30)
sigma <- sqrt(pow(invsigma2, -1))

## hyperpriors

mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))

}"





params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")



the_data <- list("y" = pt_bias_final$Bias,
                 "PartyIndex" = pt_bias_final$Partyiid,
                 "HouseIndex" = pt_bias_final$Companyiid,
                 "N" = length(pt_bias_final$Bias),
                 "J" = length(unique(pt_bias_final$Partyiid)),
                 "H" = length(unique(pt_bias_final$Companyiid))
                 )


params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")


mod = jags.model(textConnection(mod_string), data=the_data, n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 69
##    Unobserved stochastic nodes: 14
##    Total graph size: 260
## 
## Initializing model
update(mod, 1e3)

mod_sim = coda.samples(model=mod,
                       variable.names=params,
                       n.iter=5e3)

mod_csim = as.mcmc(do.call(rbind, mod_sim))


(pm_params = colMeans(mod_csim))
##   invsigma2     invtau2          mu mu_house[1] mu_house[2] mu_house[3] 
##  0.16289180  2.86648265 -0.29582529 -0.23729588 -0.36226506 -0.30358782 
## mu_house[4] mu_house[5] mu_party[1] mu_party[2] mu_party[3] mu_party[4] 
##  0.04901296 -0.89755322  0.13634614 -0.20065987 -0.14169771  0.01642649 
## mu_party[5] mu_party[6] 
## -1.03181489 -0.76678259

Model with house, party and election type

mod_string2 = " 
model {
## sampling
for (i in 1:N){
   y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] + mu_elec[ElecIndex[i]] , invsigma2)
}

## priors
for (j in 1:J){
   mu_party[j] ~ dnorm(mu, invtau2)
}

for (h in 1:H){
   mu_house[h] ~ dnorm(mu, invtau2)
}

for (r in 1:R){
   mu_elec[r] ~ dnorm(mu, invtau2)
}


invsigma2 ~ dgamma(1.5, 15)
sigma <- sqrt(pow(invsigma2, -1))

## hyperpriors

mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))

}"



set.seed(115)


params = c("mu_party", "mu_house","mu_elec", "mu", "invtau2", "invsigma2")

the_data2 <- list("y" = pt_bias_final$Bias,
                 "PartyIndex" = pt_bias_final$Partyiid,
                 "HouseIndex" = pt_bias_final$Companyiid,
                 "ElecIndex" = pt_bias_final$Eleciid,
                 "N" = length(pt_bias_final$Bias),
                 "J" = length(unique(pt_bias_final$Partyiid)),
                 "H" = length(unique(pt_bias_final$Companyiid)),
                 "R" = length(unique(pt_bias_final$Eleciid))
)

mod2 = jags.model(textConnection(mod_string2), data=the_data2, n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 69
##    Unobserved stochastic nodes: 18
##    Total graph size: 354
## 
## Initializing model
update(mod2, 1e3)

mod_sim2 = coda.samples(model=mod2,
                       variable.names=params,
                       n.iter=5e3)

mod_csim2 = as.mcmc(do.call(rbind, mod_sim2))

traceplot(mod_sim2, col = 1:3)

gelman.diag(mod_sim2)
## Potential scale reduction factors:
## 
##             Point est. Upper C.I.
## invsigma2            1          1
## invtau2              1          1
## mu                   1          1
## mu_elec[1]           1          1
## mu_elec[2]           1          1
## mu_elec[3]           1          1
## mu_elec[4]           1          1
## mu_house[1]          1          1
## mu_house[2]          1          1
## mu_house[3]          1          1
## mu_house[4]          1          1
## mu_house[5]          1          1
## mu_party[1]          1          1
## mu_party[2]          1          1
## mu_party[3]          1          1
## mu_party[4]          1          1
## mu_party[5]          1          1
## mu_party[6]          1          1
## 
## Multivariate psrf
## 
## 1
dic = dic.samples(mod2, n.iter=1e3)


(pm_params = colMeans(mod_csim2))
##   invsigma2     invtau2          mu  mu_elec[1]  mu_elec[2]  mu_elec[3] 
##  0.20759134  2.16131742 -0.31059633 -1.38158267 -0.26186068  0.63429654 
##  mu_elec[4] mu_house[1] mu_house[2] mu_house[3] mu_house[4] mu_house[5] 
## -0.47479163 -0.21559704 -0.41189908 -0.24463000  0.04654078 -0.96681012 
## mu_party[1] mu_party[2] mu_party[3] mu_party[4] mu_party[5] mu_party[6] 
##  0.18443221 -0.33279075 -0.26750023  0.30600927 -1.21179744 -0.78663631
mod_string2 = " 
model {
## sampling
for (i in 1:N){
   y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] + mu_elec[ElecIndex[i]] , invsigma2)
}

## priors
for (j in 1:J){
   mu_party[j] ~ dnorm(mu, invtau2)
}

for (h in 1:H){
   mu_house[h] ~ dnorm(mu, invtau2)
}

for (r in 1:R){
   mu_elec[r] ~ dnorm(mu, invtau2)
}


invsigma2 ~ dgamma(1.5, 30)
sigma <- sqrt(pow(invsigma2, -1))

## hyperpriors

mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))

}"



set.seed(115)


params = c("mu_party", "mu_house","mu_elec", "mu", "invtau2", "invsigma2")

the_data2 <- list("y" = pt_bias_final$Bias,
                 "PartyIndex" = pt_bias_final$Partyiid,
                 "HouseIndex" = pt_bias_final$Companyiid,
                 "ElecIndex" = pt_bias_final$Eleciid,
                 "N" = length(pt_bias_final$Bias),
                 "J" = length(unique(pt_bias_final$Partyiid)),
                 "H" = length(unique(pt_bias_final$Companyiid)),
                 "R" = length(unique(pt_bias_final$Eleciid))
)

mod2 = jags.model(textConnection(mod_string2), data=the_data2, n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 69
##    Unobserved stochastic nodes: 18
##    Total graph size: 354
## 
## Initializing model
update(mod2, 1e3)

mod_sim2 = coda.samples(model=mod2,
                       variable.names=params,
                       n.iter=5e3)

mod_csim2 = as.mcmc(do.call(rbind, mod_sim2))

traceplot(mod_sim2, col = 1:3)

gelman.diag(mod_sim2)
## Potential scale reduction factors:
## 
##             Point est. Upper C.I.
## invsigma2            1          1
## invtau2              1          1
## mu                   1          1
## mu_elec[1]           1          1
## mu_elec[2]           1          1
## mu_elec[3]           1          1
## mu_elec[4]           1          1
## mu_house[1]          1          1
## mu_house[2]          1          1
## mu_house[3]          1          1
## mu_house[4]          1          1
## mu_house[5]          1          1
## mu_party[1]          1          1
## mu_party[2]          1          1
## mu_party[3]          1          1
## mu_party[4]          1          1
## mu_party[5]          1          1
## mu_party[6]          1          1
## 
## Multivariate psrf
## 
## 1
dic = dic.samples(mod2, n.iter=1e3)


(pm_params = colMeans(mod_csim2))
##   invsigma2     invtau2          mu  mu_elec[1]  mu_elec[2]  mu_elec[3] 
##  0.18720907  2.30782506 -0.30919797 -1.26847525 -0.27761760  0.57370895 
##  mu_elec[4] mu_house[1] mu_house[2] mu_house[3] mu_house[4] mu_house[5] 
## -0.48394649 -0.20971232 -0.39773000 -0.24919981  0.03576236 -0.93825866 
## mu_party[1] mu_party[2] mu_party[3] mu_party[4] mu_party[5] mu_party[6] 
##  0.16292985 -0.32405373 -0.25471263  0.25056884 -1.15488402 -0.76468216